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A novel approach to the problem of the slit correction in small-angle x-ray scattering is presented, 
based on a matrix inversion method. The integral equation for the slit correction is written as a Volterra 
equation of the first kind. This equation is reduced to a system of simultaneous equations, expressed 
in matrix form. The order of the matrix is equal to the number of experimentally determined points. 
To obtain these equations, one has to expand the unknown functions in Taylor series around each 
point to be subsequently determined. There is, however, a difficulty inherent in this method due to the 
fact that most of the series expansions of the function to be determined lead to a strong numerical 
instability. However, a general method is developed, which enables us to find shifting operators leading 
to numerically stable systems of equations. The "unsmearing" of the experimental data is then performed 
by standard matrix inversion procedures. 
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1. Introduction 

In a previous paper [1] ' methods for solving the integral equation for slit correction in small 
angle x-ray scattering experiments were proposed. This equation relates the experimentally 
determined intensity of scattered rays to the true scattering intensity as follows: 



7(|) r~ y /(y)r(v^i)<fr (1) 

Jx Vy 2 — x 2 

In eq (1), I(x) is the experimentally determined intensity of the scattered light, I(x) is the true 
scattering intensity, x is proportional to the scattering angle [1], and W(x) is the experimentally 
determined slit weighting function. The dependence of W on the slit collimation system is widely 
discussed in references [2-5]. It may be given either in analytic form or only for discrete values of 
its argument. 

I(x) is to be calculated from eq (1), which is a Volterra integral equation of the first kind. 
In Reference [1] explicit and implicit solutions have been proposed. The explicit solution is based 
on transforming eq (1) into a Volterra equation of the second kind, which has for its inhomogeneous 
part a constant term, and whose kernel involves an integrated expression of the derivative of the 
experimentally determined slit weighting function W. This solution requires numerical differentia- 
tion of the experimental data, a process which may introduce serious errors. Differentiation is not 
required in the implicit solution method. Rather, the experimental data enter into the solution 
directly. 
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This paper will deal only with the implicit solution method. In that case, eq (1) is first reduced 
to a system of simultaneous linear equations. These equations are then solved by inversion of the 
resulting triangular matrix. The solution leads to calculated approximate values of /(:*;), given on a 
set of TV discrete points x=0, x = Ax, . . ., x= (N-l)Ax. These are the same points for which 
values of I(x) must be given by experiment. 

One of the difficulties in applying this procedure directly is that it can lead to numerical 
instability. However, one can find an interpolating polynomial for I(x) such that a system of 
simultaneous equations is obtained which is numerically stable. In this paper will explore ways 
to find an accurate interpolating polynomial for the purpose of our computation. However, the 
problem of finding a suitable polynomial can be altogether avoided if another, indirect method is 
employed. That method consists in first obtaining a zeroth-order approximation to the solution, 
based on a single-step method and which is known to be always numerically stable. Subsequent 
improvements to the solution are then computed, using standard perturbation methods to obtain 
subsequent corrections to the zeroth-order approximation. 

The first of these two approaches involves application of multi-step methods which are based 
on existence of polynomial which interpolates the function I(x) at discrete points on a given 
interval of x. The second approach is based on a single-step method, which serves as a starting 
point for the perturbation method. Both approaches are included in a computer program that 
solves the problem of slit correction in small angle x-ray scattering. Four different choices for 
the weighting function (W in eq (1)), have been used in that program: 

1. Constant weighting function (infinite slit) 

2. Gaussian 

3. Trapezoid, and 

4. Data weighting function. 

The details of these computations and an analysis of the results will appear in an NBS Technical 
Note, to be published. In this Note, the Fortran program that performs the required calculations will 
be described in detail. 

2. Algebraic Formulation of the Problem 

The kernel of the integral eq (1) has a singularity aly = x. In order to deal with this singularity, 
we follow the procedure employed in solving Abel's integral equation. Both sides of eq (1) are 
multiplied by xdx/Vx 2 — u 2 and integrated over x from u to infinity. The transformed equation,* 
after exchanging the orders of integration is 



I I(x)G(u, x)dx= I I(x)K(u, x)dx 

Ju Ju 



(2) 



K(u,x) = 2x [ V "- , " T7 = (2a) 



ith 

'v^^ W(t)dt 

V* 2 - u 2 - t 2 
and 

G(u, x) = x\ Vx 2 — u 2 

The total number of experimental data for I(x) is N, which are given at equally spaced intervals Ax. 
The range of x is determined in such a way that putting I(x) equal to zero for all:* ^ NAx does not 
affect the solution. In the physical problem we are investigating, the assumption that I(x) = for 
x ^ NAx implies that I(x) = for the same values of x. Therefore, both the upper limits of integrals 
in (2) are replaced by NAx. Thus, the function I(x) represents a discrete set of function values /*, 
i= 1,2, . . ., N. Both I(x) and I(x) may therefore be viewed as TV- dimensional vectors. The first 
components Fi and h of these vectors correspond to a = 0; i = N corresponds to u = (N — 1 ) Ax. 
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It should be remarked, that consistent with the above restrictions on the functions I(x) and 
I(x), the method of solution presented in this paper is not necessarily restricted to the particular 
representations of K(u, x) and G(u, x) , which in the problem being investigated are given by 
eqs (2a) and (2b). Hopefully, the method presented here is general and might be applied to other 
physical problems which can be simulated by eq (2). 

The left-hand side of eq (2) is calculated from known quantities. The results are represented 
by the vector h with components 



[NAx m 

l i;= I{x)G(u, x)dx. 

J(/-!)A.r 



(3) 



In our particular problem, G(u, x) is singular at x = u and is given by eq (2b). For this reason, 
a direct numerical evaluation of the integral (3) will not be practical, since we will have to sub- 
divide the integrand very closely at its lower end in order to achieve sufficient accuracy. To evaluate 
(3) we replace the integral by summation as follows: 



/, 



x xl(x)dx _^ P Ar xl(x)dx 

u y/ x --u- ~ pi J o- 1 >A* y/x' 1 - (i-l)-Ax' 2 ' 



Next,/(#) is expanded in a Taylor series around x=jAx: 

xdx 



N fjXr r 

tl Jo-da* L 



I(jbx) + (x-jkx)Ai+- (x-jAx)-A 2 + . 



Vx--(i-])' 2 Ax- 



(5) 



A], A 2 , etc., are the first, second, etc. derivatives of I(x) evaluated at x=jAx. These derivatives 
can be approximated by one of several difference relations. We may use the following finite dif- 
ference expression: 

Ai=(F j+l -lj)IAx, A. 2 =(I j+2 -2I j -i + I j )l(Axy. 

Only two derivatives are used in the Taylor series in (5); I N+i and I N+2 are set equal to zero. After 
performing the indicated integrations and rearranging eq (5), the result is written in the form 

hi=^yt,jlf, lj = (orj>N. (6) 

or h= AI. A is an upper triangular matrix with entries the coefficients yij. The matrix elements 
ytj, have different representations depending on the following three cases: (I) j—i = Q, (2) j—i 
even, and (3) j—i odd. In addition, in order to close the system, one has to consider, separately, 
the cases j=N and j = N—l. This is because the last integral, if taken from x = (iV — 1) Ajc, cannot 
include second-order differences if forward extrapolation is employed. Also for the case i=l(u = Q) 
we have to avoid the expression (i — l)log(i — 1). Details will be given in an NBS Technical Note 
describing the computer program. 

We now expand the r.h.s. of the integral equation (2) in terms of the discrete values /,, as 



/; 



v 



I(x)K(u,x)dx=Yb iJ I j ; lj = for j>N. (7) 

-1)Ajt ^ 



J=i 



to obtain a set of simultaneous equations 



Yhth-hi (8) 
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bu j(i=l, 2, . . ., N; j=i, i+1, • . ., AQ are coefficients which we have not yet determined. 
These coefficients depend on the functional form of the kernel K(u, x), and, therefore, on the as- 
sumed form for the weighting function W of eq (1). In addition, they also depend on the coefficients 
of the interpolating polynomial which approximates the function I(x). 
The solution of eq (8) can be expressed in matrix form: 

I = B h. (9) 

I is the TV-dimensional solution vector, h is an /V-dimensional vector whose elements were already 
computed from the experimental data and the indicated numerical integrations. The inversion 
of an upper triangular matrix B can be rapidly performed by a back substitution [5], starting with 
the last equation. With this method, the elements cij of the inverse matrix B _1 are obtained by 
first calculating the elements d,j from the following recurrence equations: 

and then a, j = dj, j/bi, j . 

For certain forms of K(u, x) , the kernel of the r.h.s. of eq (2), there is a repetitious pattern in 
the representations of the elements djj. In this case, the number of different elements in the matrix 
B is less than (1/2)N(N + 1). In particular, if K{u, x) is independent of u (e.g., when W of eq (2a) 
is constant and therefore, K{u, x)=2x) the elements of the matrix B (after multiplication by the 
diagonal matrix 

depend only on \j — i\. This has an important practical aspect: computer computations with large 
matrices require storage capacity which is often not available. It is then necessary to employ vectors 
bj (formed separately for each row of the matrix B), or, to compute bij each time it appears in 
computations, separately. The problem is, that these computations will be repeated many times. 
However, for this case, only the first row of the matrix, in form of a vector b \ needs to be computed, 
and the problem of storing large matrices is avoided. 

Any algorithm in which the calculated value of the increment I\—h + \ depends only on/ i+ i and 
on the step Ax, is called a one-step method. In other words, an expression for It may be written 
in a closed form involving explicitly only 7,+i Of course, U+\ depends, in its turn, on 7/ +2 , etc. In 
a multi-step method of multiplicity /, the increment If — Ii+i depends also on /j+2, • • .,/?+/• It 
therefore follows, that whenever one approximates I(x) in eq (1) by an interpolating polynomial of 
order higher than one, the method of calculating U in terms of already computed /,-+i, /j+2, etc., 
will be defined as a multi-step method of multiplicity which is equal to the order of the interpolat- 
ing polynomial. In a multi-step method, a special starting procedure is required to find the "initial" 
or "starting" values of In, In-i, • • *,In-i- We calculate these starting values by trapezoid rule as 
follows: 

1 C(N-l)Ax 1 CNkx 

hs-i=j: (lN+lN-i)K[(N-2)bx 9 x}dx+±I N K[(N-2)Ax,x]dx 

2 J(N-2)Ax 2 J(.V-i)Aj» 

N _ l = I x _ l \ K[(N-l-l)bx,x]dx + I N -M K[(N-l+l)tkx 9 x]dx 

J(X-I-\)Ax J(N-l-l)Ax 

CNAx 

. . . + I N K[(N-l-l)Ax,x]( 

J(N-2)Ax 
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2h 



]dx 



We have / equations and /+ 1 unknowns. However, a physical assumption was made that setting 
h = for i > N, implies that /* = for i > N. For this reason, 

I N =2h N I \ J K[(N-\)Ax,x]dx. 

I J (N-l)bx 

Let us divide the integral h, 

hi=\ I(x)K[(i-l)Ax,x]dx 

J (/-1)A.r 

into smaller integrals, each with a range of integration equal to lAx (except for the last integral). 
For this purpose, let us define 

Kij=\ I(x)K[(i-l)Ax,x]dx V=Ui + h- • .,A/-1) (11) 

J U-DAjc 

so that Ki, i = hi. 

hi is the sum of differences of the Kij as follows: 

hi= (Ki,i — Ki, i+l) + (Ki,i+l — Ki,i+2i) + - • •• 

(In this expression, the number of bracketed terms is equal to the nearest integer not exceeding 
(N — i)IL In addition, a term is added to close the system so that hj = K,j.) Next, in each of the 
difference equations 



nj+i-D&x 
Ki,j — Ki,j+i= I I(x)K[(i — 1)Ajc, x]dx 

J(j-DAx 



(12) 



we replace I(x) by an interpolating polynomial having the values Ik = I(xk) on a set of points xi<. 
These points might be taken as Xk= (/— 1)A#, . . ., (j — I— I) Ax. In forward interpolation, the 
interpolating points are the p points j—\ + m, j+rn, . . ., j+ m + p — 1, with p^l. m denotes 
the point from which the differences are computed. We will always have m < /, so that the point 
of reference lies within the range of the integral (12). The polynomial Q(x) which approximates 
and replaces I {x) in eq (12), in terms of forward differences [6], is 

(?(*)= J HxlAx)-j-m+l\teI j+m . (13) 

«-° \ q I 

In eq (13), A is a forward difference operator. Thus, for #=1, A// +m +i = / +/w+1 — //+ m , etc. 
Substituting eq (13) in eq (12), integrating and rearranging, the result is formally written as 

K iJ~ Ki,j + l = (*ijIj + m + (XiJ+\Ij + m+l+ • • • + «i, j+plj+m + p (14) 

«i,j+ft(A; = 0, 1, . . . , p) are the coefficients resulting from the above steps taken in order to 
approximate eq (12) by eq (14). 

In the formal solution, written in matrix form as 1= B _1 h, the elements of the matrix B can 
be expressed in terms of aij of expansion (14). 

Before proceeding, let us stop and see what will happen to our computations with I(x) of eq 
(12) being replaced by polynomial Q(x), for the specific case of m = and l = p = 2. In this case, 
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all three interpolated values of U are within the range of the integral of eq (12). Let us assume that 
K(u, x) oc x. This is the case when the weighting function W is constant (W= constant is descrip- 
tive of an infinitely-wide slit in small angle x-ray scattering). For i > 1, the elements of a typical 
row of the B matrix are found to be 

i 4 2 

bi,i = - ; bi,i+i=- (i + 1); bi,i +2 = - (i-f 2). 



The coefficients aij of eq (13) are: 



4 1 

3— 3 U + 1);a = 3 



a/,/+i=- (i + 1); a=- (£ + 2). 



The absolute values of the elements of the inverse of the B matrix can be shown to grow exponen- 
tially as one moves away from the diagonal elements. Consequently, the numerical calculations 
are affected by strong numerical instability. 

3. Stability Conditions for Integral Equations of First Kind 

3.1. General Theory 

Consider first the difference eq (14) written for k/, / — Ki,i+i and for p=l. The associated 
homogeneous equation, obtained by setting 

Ki,i — Ki, ;+/ = 0, is 

«/, /// + //* + «/,,■ + ! Ii+m + i + • • .+OLi,lIi+m+l =z 0. (15) 

Solution of (15) can be found for every choice of starting values 7 \ , /\-i, . . ., /.v-z+i. (These are 
computed using the trapezoid rule and the one-step method.) In particular, if these starting values 
are all equal to zero, then there corresponds a trivial solution /, = for every i. 

Consider first a simplified homogeneous equation for the starting values of I, with constant 
coefficients, replacing eq (15), 

ao/.v-/+ai/.v-/+i + • • • + O!//a™0. (16) 

We determine the conditions for convergence of solution of eq (16) for /V— > oo, while retaining 
finite value for /VAx. That is, the intervals between successive points becomes smaller and smaller 
as TV increases, and, as N-*™, Ax— »0. If the multi-step method for solving eq (16) converges, 
then it will converge also for the special value of Is — 0. 
Therefore, 

lim /„ = (), il = N, /V-l,. . ., N-l (17) 

is a condition for convergence of multi-step method for solving eq (16). Accordingly, if condition 
(17) is not satisfied, the solution of the general homogeneous eq (15) will diverge in the limit of 
Ax— »(). This divergence implies the divergence of the solution of the inhomogeneous eq (14). 
The procedure to determine the conditions for convergence of solution of eq (15) in the limit of 
Ax^ °° (or A^-^ oo) , is identical with the one employed in connection with numerical solution of 
difference-differential equations [7], and proceeds as follows: The homogeneous eq (15), with con- 
stant coefficients, is expressed as: 

at/i+m-f-ai+i/i+m+i+ • • • + ai+Ji+ m +t=0. (18) 
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In eq (18), we substitute 

Equation (18) becomes after division by common factor r)~ (i+m+t) : 

Therefore, r\ must be a root of the characteristic polynomial ^(17): 

0(T/) = a,y + a/ +1 i7 / - 1 + . . . + o /+ , = 0. 



(19) 



(20) 



We require that the multi-step method be convergent. The necessary condition for the con- 
vergence of the general multi-step method for our integral equation of the first kind is the same 
as the condition for the convergence of difference equations of first order [8], and is stated by the 
following: 

THEOREM: A necessary condition for the convergence of the multi-step solution (14) for in- 
tegral equation (1) is that the moduli of all roots of the associated polynomial '(8(17) do not exceed 
unity, and that the root of modulus 1 is simple. 

The condition imposed on 0(17) is also the condition of stability [9], Violation of this condition 
leads to "strong" numerical instability. This instability manifests itself by the error growth at a 
fixed point x = (i — 1)A* as the spacing between the points decreases. Strong instability is therefore 
characterized as a pointwise phenomenon, as it involves error growth at fixed points as the interval 
between them decreases. In contradistinction to this, a "weak" instability, to be discussed sub- 
sequently in connection with even-odd oscillations in computed values of//, is associated with the 
truncation error, and exists even when the method employed is convergent. Since the weak insta- 
bility is associated with the total range of x, while keeping the interval between successive points 
constant, this instability is also known as a stepwise phenomenon [7, 10]. 

Let us demonstrate the concept of strong instability by a concrete example. Assume that 
/ = 2 and m = 0. Also, let K(u, x)ccx, so that Kij — Kij+2 = /?, — /?;+ 2 . We substitute eq (13) into eq (12) 
and integrate. Since the r.h.s. of eq (12) is a function of j only, ctij, the elements of the matrix B, 
are shown to be equal to ja\j-t\ for j ^ i and to zero for j < i. ol\j-%\ is a function of/'—/ only. If, in 
addition, p=/, then, for i> 1, 



hi — h, +2 = aj; -h ai/,4 1 + otzl 1+2 



(21) 



and 



00 = 0,, i 



'- «/, 1 



; lOCo, (X] —OLi, j+1 



ia-i. 



l(X\ 



For a particular case when 0:0 = 0:2 (as in the example discussed at the end of chap. 2), the ele- 
ments b*i,j of the inverse matrix B ' can be expressed directly in terms of the roots of the charac- 
teristic polynomial ty{j)), eq (19). One obtains 



ib* Uj = 
U>i) 



and b*i,i= 1//. 



ZOo 



-h Vor — 4oj 



2oo 



- v ot\ — 4oiOo I 



(22) 
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In eq (22), k=j—i. When a 2 9^ «o, the relationships between 6*/j and the roots of eq (19) are 
more involved. However, for sufficiently large A, the leading term in the expression for b*i,k+i can 
be shown to be equal to 



-OL\ 



2a 
For the case of K(u, x)ccx, the following equation is derived 

/»-/i + «=2 &* i , i [Aj-(i/(i + 2))A j+2 ]. (23) 

When k becomes infinitely large U will remain finite only if lim b*ij converges. This happens 

only when both expressions in parentheses in eq (22) do not exceed 1 in their absolute values. This 
implies that the moduli of the roots of the equation 

0(-n) = r]--\- aiT)-\- a* — 

do not exceed unity. This simple example demonstrates the intimate connection between the root 
condition and the stability of the multi-step method. 

In the general theory of stability of the difference differential equations, the necessary and 
sufficient conditions for the convergence of the multi-step method require, in addition to the root 
condition, also the satisfaction of the consistency condition [9]. This condition requires that, in 
order that an expansion in discrete variables (e.g., our expansion eq (14)) be a good approximation 
to the original equation which it simulates, certain restrictions must be imposed on the coefficients 
in such expansion. In our problem the consistency condition requires that the sum of the coefficients 
in eq (18) be equal to /: 

r=i+l 

2 « r =Z. (23) 

r=l 

Since this condition is always satisfied, the root condition becomes both the necessary and the suf- 
ficient condition for the convergence of the multi-step method. We also see why the method for 
the determination of the starting values of //, based on a one-step method, is always stable and 
convergent: The polynomial 0(r)) is then of the first order with its root equal to 1. 

3.2. Determination of Stability Conditions 

Our purpose is to obtain an approximation to the integral 

I(x)K[(i-l)Axjc]dx = Ki,i-Ku+i (24) 

J(i-l)Ax 

in terms of finite differences, which satisfies the stability conditions. These differences can be re- 
ferred to the point x = (i — 1 ) Ax , to the point x = iAx , and so on. Let the point of reference be located 
at x= (i+ m— 1)A#, /m = 0, . . ., /. We define a shift (displacement) operator E by 

EWj-i. (25) 

Applying the shifting operator to the point of reference, one obtains 

I(x) = E- t (I i+m ). (26) 

t is a new variable, defined by t = xjAx — (i -+■ m — 1 ) , and it measures the distance from the point of 
reference. Therefore, 
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f(l-m) 

Ku-Ki, i+ i = Ax \ K(i*,i*+t*+m*)E-*lI t+m - 1 ]dt. (27) 

J—m 

In eq (27) and in subsequent equations, i*=(i — 1)A#, /*=jAx, and m*=mAx and /v(i*, /'* + 
t*-\- m*) = K[(i— 1)A.%, (i + 1 + m — 1 )Ax] . The operator E can be expressed in terms of one of 
the three difference operators, A, V, and 8 which are defined by 

Mj = lj+\ — Ij (forward difference operator), 

AIj= Ij — Ij-i (backward difference operator), 

8Ij = /j+i/2 — //-i/2 (central difference operator). 

These operators are related to E through 

A=l-E, A=E- 1 -l,8=E 1 / 2 -E- 1 / 2 . 

For the backward difference operator A we have the following eq: 

fl-m 

Ki, , - K U i+i= Ax K(i* 9 1* + t*+ m*) [( 1 - A)-'/i + «-i]cfc. (28) 

For forward difference operator A we have 



ri-m 
K U i-Ki, i+ i = Ax K(i* 9 i*+t*+m*) [(l + A)7 <+w -i]&. 

J— Wl 



(29) 



In order to derive an equation for k,-, ,■ — k;, ;+/ for central differences in terms of difference 
operator 6, we have for express E in terms of d. Let E=e D , and, therefore, 5=2 sinh (D/2). Thus, 

D = 2sinh-> (8/2) = 8-|^ + H^+. . . • 

E "'= 1 - / ( 8 -i^ + ---) + f( s -3i + ---) 2 - ( 3 °) 

Retaining terms involving no higher operators than 5 2 , we have E - ' = 1 — dt + t 2 8 2 /2. 
Therefore, 



f 1 -™ • r/ / 2 /5 2 \ 

ici,i-Ki > < +2 =A« £(**, i* + **+m*) n-8t+-f\I i+ 



dt • (31) 



Let us now investigate a few cases with simplified kernels, and find out which ones lead to numer- 
ically stable solutions. For the kernel K(i*, i* + 1* + m*) we again adopt the form which is descrip- 
tive of the constant weighting function in small-angle x-ray scattering, that is, K(i*, i* + t* + m*) 
= (i-f t-\- m— l)Ax. For the first case let us take the form is which the approximating polynomial 
agrees with I(x) at x= (i — 1)Ajt, . . . , (i+l— I) Ax. We take 1 = 2, and m = 2. In order that all 
points will fall within the range of the integral, we must use the V (backward) operator. Consider 
i> \. The root-determining equation is 1/3 (x 2 + 4jt -f 1) = 0. (The coefficient of x 2 is the coefficient 
of //, the coefficient of x is the coefficient of /,-+i, and the last term is the coefficient of /j+2.) We 
detect immediately that this form leads to a divergent solution. (Root condition is not satisfied.) 
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Let us now consider the case with central differences centered about the mid-point of the in- 
tegrand. Let 1 = 2, and m=\. From eq (31), the following root-determining equation is derived: 



1/3*2 + 4/3*+ 1/3 = 



(32) 



and the method is unstable. It is interesting to notice that this particular case, which is the most 
obvious one since it is symmetric and does not involve odd powers of t (in the integrand of eq (31), 
leads to an equation which is equivalent to the extended Simpson rule for numerical integration. 
The same result would be obtained if we evaluated the integral as 

CNkx 

f(iAx,x)dx = ll3(fi, i +if i , i+1 + 2f iM2 + . . .)Ax (33) 



f 



with/(iAx,*) = K(i\x,x)I(x) , etc. 

At this point, we are led to the conjecture that the most straightforward cases, namely the ones 
that involve points located within the range of the integral, lead to unstable solutions. If the above 
restriction is removed, stability conditions may be satisfied. For example, in the case with m = 2, 
(the finite differences refer to the point x= (i + l)Ax) we have, for / = 2 and the A operator, two fixed 
points which we located outside the range of the integral. The root condition can be shown to be 
satisfied. 

The method that we find to be the best for our purposes is the one in which the A operator is 
used and the absolute values of the coefficients of A', r= 1,2, . . . ,p, are made as small as possible. 
This can be achieved by selecting the point of reference to be located in the middle of the inte- 
grand, i.e., by taking m = 1/2. In this case, one is justified in retaining only the first few terms in the 
expansion of (27). We have, for 1 = 2, 



Ki,i+2 Z 



I 



Ax\ K(i*,i*+t*) 



1-;A+^ (f-l)A 2 -. 



hdt. 



(34) 



For K(i*,i* + t*) = (i+t)kx, and i > 1, the coefficients of h+j(j= 0, 1 or 2) in the expansion (34), 
after integration, are identical to the ones employed in Nystrom's method [11]. Since the root de- 
termining equation is 



7/3x 2 -2/3x + 1/3 = 0, 



(35) 



the root conditions are satisfied. We found that, of all the multi-step methods, the above method is 
the most accurate for our purposes. Part of the computer program for the slit correction in small- 
angle x-ray scattering is based on this algorithm. 

In conclusion, to insure stability, the differences in the expansion of the shifting operator should 
be referred to the point (/*•) to be calculated and forward differences should be used. We also con- 
jecture that the central difference method, which is the most accurate one in solving ordinary 
difference differential equation, is always numerically unstable. 

The following qualitative considerations indicate that the stability conditions are only slightly 
altered by the actual form of K[(i — i) Ax, x] and that a non-constant weighting functio n relaxe s 
the condition for the attainment of the stability conditions. Consider K{u, x)=xF{ V# 2 — u 2 ) , 
(which is the form used in actual computations of slit corrections). F is related to the weighting 
function W by 



F(Vx~ 2 



l 



\Zx 2 -U* 



W(y)dy 



Vx^- 



■f 



(36) 
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Now, let us form a matrix of bt, j, where btj are the coefficients in the approximation of 



rxAx 

Ju 



xI(x)F( Vx 2 — u 2 )dx by an expansion V bijlj. 

Consider, for example, a trapezoid form of W(x). It can be shown that in each row of the matrix 
B= (bjj), bij will decrease as j increases, as compared with the case of constant W. This will 
result in a decrease of the ratio of the first to the last coefficients in the root-determining equation, 
and, therefore, in a decrease in the value of the modulus of the largest root of the polynomial, eq 
(20). 

Thus, the stability conditions severely restrict us in selecting the best procedure for numerical 
computations of a system of simultaneous equations which replace the original integral equation 
of the first kind. This stability problem can be almost completely eliminated if the original integral 
equation of the first kind is transformed into an equation of the second kind, by differentiation. 
Then eq (21) becomes 

-K{u,u)I(u)+ I(x)^K(u,x)dx=^< (37) 

J,, du du 

\{K{u, u) = 0, the differentiation is repeated, until one obtained 



d" 
^K(u.x) 



#0 
x = 



When eq (37) is replaced by a set of simultaneous linear equations, the presence of the first term 
in the l.h.s. of eq (37) drastically relaxes the restrictions imposed by stability consideration: 
The reason is that the presence of this term increases the values of the diagonal terms in the B 
matrix. This, in turn, tends to decrease the modulus of the largest root of eq (19). For this reason, 
it is preferable to transform eq (2) into eq (37) whenever / (u) and K(u, x) are given as analytical 
functions. However, in our problem these functions are given as sets of points based on experi- 
mentally determined data. Their differentiations would introduce a further uncertainty in the 
computed results. The method of solution we developed is to be preferred only if the numerical 
differentiations of experimental data are to be avoided. 

4. Perturbation Methods 

One would expect that, of all difference approximations to derivatives, the ones based on cen- 
tral differences would lead to the most accurate results. This would be the conclusion based on 
formal identities between the derivative operator D r , defined by 

Ul \ dx<- )x=(i-l)Ax (38) 

and the various shifting operators A r , V r , and S r . 

From E = e D (in units employed in this work, based on arbitrary setting of Ajc = 1) we have, to the 
first approximation, 

D = (log E)' = A^l -'- A)s V '(l + r/2V )■ S'(l -~^j ■ (39) 

Equation (39) shows the magnitude of the error involved in replacing an r-order derivative by various 
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r-order differences. It is the method, based on central differences, which leads to the standard inte- 
gration formulas. For example, with 1 = 2 and p = 2, the coefficients in expansion (31) are the same 
as in Simpson's rule for numerical integration. With l = p=6, Weddle's rule is obtained. Unfortu- 
nately, the direct solution of the integral equation, based on the central-difference method, was 
shown to be numerically unstable. If, however, we wish to use the central difference method, we 
will have to employ other, indirect approaches to the algebraic solution of the integral equation. 
These approaches will be discussed now. 
We start with 

i=N j=N f(j-l+//2)Aj- 

2 bt,jlj= 2 h\ K[(i-l)Ax,x]dx 

j=U M, . . . j=i, i+l, 1+21, . . . J(j-l-ll2)*x 

+ 1/2[/j+i-2/j + //-i] I 2 X (x-j-l) 2 K[(i-l)Ax, x]dx+ . . .=ht • (40) 

In matrix form, BI = h. Expand B = B () + B2 + B 4 + . . . . The elements of the matrix Bo, bo(i,j), 
are 

f(j-\ + U2)Ax 

bo(iJ)= K[(i-l)bx,x]dx 

J(j-l-ll2)Ax 

for 7= i, i-h/, i + 2/, . . . and 6o(/, j) =0 for all other values of 7. B 2 is the matrix of 62(1,7), and 

f(j + l/2-l)Ax 

b 2 (ij)=- [x - (j ~ I) hx] 2 K[(i- I) &x,x]dx for j=i, i+l, i + 2l , 

J0'-//2-1)Aj- 

r (J-i+II2)Ajc 
b 2 (ij)=+ (x-jAx) 2 K[(i-l)Ax, x]dx 

JU-\-U2)toc 



0-1+ 1/2) &r 
J(j.-l-lj2)kx 

for 7=1+1, j + Z+1 . . . , 



b 2 (i,j) = H2 f° l + ll2)1X [x-(j-I)hx] 2 K[(i-l)Ax,x]dx 

J(j-l-l/2)lx 

for j=i— 1 and 62(1, j) = for all other values of 7. Similarly expressions for 64(^,7), etc., can 
be formulated. Let the zeroth-order solution be 

Bol (0) = h. (41) 

Application of perturbation theory leads to the following: 

B = B» + XB 2 + X 2 B 4 + . . . (42) 

I = I<"> + \I< I > + A 2 I< 2 >+. . . . (43) 

Equating equal powers of X in equation BI = h leads to the following expressions for I (1) and I ( ->: 

/(i) = _B-iB 2 I (0) (44) 

and 

P> = -B -' [B 4 I (0 > + B 2 I (1) ]. (45) 
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Since the zeroth order solution by definition, satisfies the conditions for numerical stability, the 
elements of the matrix B^ 1 are bounded. Specifically, for K[(i — 1) Ajc, x]=x and / — 2, B^ 1 is 
is of the form 



B, 



the matrix B 2 is 



1/2 





-1/2 











1/4 





-1/4 











1/6 





-1/6 



B 2 = l/6 



2 


4 


-6 


8 


-10 


2 


-4 


6 


-8 


10 





3 


-6 


8 


-10 








4 


-8 


10 











5 


-10 



We obtain for the correction term, I ( . l \ the following relation: 



/ci) = _l/12 [/Co) r 



.2/(o) + /(o) 1 
/ /+i J 



(46) 



Thus, the correction to the zeroth order approximation is evaluated by calculating the second order 
differences, 5 2 / (0) of the zeroth-order solution, and subtracting 1/12 8 2 I {0) from/ (0) . 

The above scheme does not exhaust all possibilities in applying the perturbation method. 
One can use the zeroth-order solution based on a trapezoid rule. The B 2 matrix will then be con- 
structed differently than the one given in the above example. However, a typical row of the inverse 
matrix B 1 (again, calculated for K[(i — l)Ax, x] = x and for i> 1) will be. 

-[1,-2,2,-2,. . .] 



resulting in even-odd oscillations in the computed values of//. This is also known as weak numerical 
instability. This weak instability can be greatly reduced, if certain averaging procedures are 
adopted. For example, one can replace, / i (1) = V /3*(i, j)Ij, where f3*(ij) are the elements of ma- 

_ J' 1 

trix B^ 1 B 2 , by I (1) =DI, where 8(i, j) , the (i, j) element of the matrix D is given by d(i, j) = 



\ [/»•(«. 



j)+p*(i+l,j) 



(47) 



In summary, two methods are at our disposition: 

(1) A multi-step method with forward differences, referred to the point being computed and 
extended to the points already computed. This method requires precalculated "starting" values. 
This number of starting points is usually equal to the multiplicity of the method. These starting 
values are most conveniently computed using a single-step method. 

(2) A single-step method used as zeroth-order approximation. Higher-order approximations 
are then computed by a perturbation method. The advantage of this procedure is in employment of 
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the more accurate central differences replacing the various derivatives. This method, if employed 
in combination with certain averaging procedures, tends to remove some of the weak instabilities 
(in form of even-odd oscillations) from the zeroth-order solution. 

5. Some Computational Details 

The two approaches outlined above, one based on a multi-step method and the other based 
on a single-step method serving as a zeroth order approximation in perturbation method, were 
employed in writing a general computer program (in Fortran) for the problem of numerical solution 
of eq (1). The first approach employs eq (34) as a starting equation for the actual computations. 
This equation is employed to the approximation that only terms up to A 2 are retained in the ex- 
pansion of the bracketed terms of eq (34). In the second approach the corrections to the zeroth 
order solution (which is given by eq (41)), were computed, by using eq (46) and also by performing 
averaging indicated by eq (47). 

In actual computations, four different forms for the slit weighting function were used: 

1. Constant weighting function (infinite slit) 

2. Gaussian 

3. Trapezoid, and 

4. Data weighting function. 

For the first three functions, parameters which describe them have to be supplied in the form of an 
experimentally determined set of values taken at a given constant interval. The kernel of the 
integral equation (given by eq (2a)) for these first three functions, is evaluated analytically. More- 
over, all the required moments of the kernel (that is, integrals of the form fx*K [(i— 1)A.%, x]dx, 
are also evaluated analytically. For the last assumed form for the weighting function, however, 
the computations of the kernel of the integral equation and its moments, have to be performed 
numerically. The time required to perform the entire "unsmearing" of the data is considerably 
longer for the case of the data-weighting function, than for either one of the other three functional 
forms of this function. 

To our knowledge, the slit-weighting functions (1— 3) are the only functions, which describe 
the actually measured weighting functions reasonably well, and which are also analytically 
tractable. 

The results of the computations were subsequently tested by an inversion of the calculation 

procedures as follows: Using the calculated values of I, the values of I were recalculated and 
compared with their original values. The accuracy of the numerical procedure was tested both 
with various trial functions for which analytical methods for solution of the original integral equa- 
tions are available, and with experimental data taken from various measurements. Both approaches 
to the solution of the starting integral equation were tried with satisfactory results. 



The author thanks R. J. Rubin for helpful discussions and suggestions in preparation of this 
publication. 
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